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1. Introduction 

A system connected to two particle reservoirs at different densities relaxes to a 
noneqnilibrinm steady state (NSS), with a particle current flowing through it. The 
description of the fluctuations of this current has recently received much attention 
[H El El @ El @ EM El HU m El EU HU C 32 ] . In equilibrium, thermodynamic potentials 
are related to exponentially unlikely fluctuations away from the average |T6] , as was 
already discussed by Einstein dt Analogously, one can construct nonequilibrium 
thermodynamic potentials from the study of exponentially unlikely current and density 
fluctuations away from the NSS ra- A theoretical framework for this approach is 
provided by the macroscopic fluctuation theory (MFT) [HI {20] [21J [22l [23] . 

Using the MFT, Akkermans and co-workers studied current fluctuations in diffusive 
systems connected to two reservoirs m- They showed analytically that for a system of 
arbitrary (but fixed) dimension, the ratio of the cumulants of the current distribution is 
independent of the shape of the system and the shape of the contacts with the reservoirs. 
This derivation is valid if both the system and the contacts with the reservoirs are 
macroscopic in size. The analytical prediction was tested by numerically calculating 
the ratio of the first two cumulants, called the Fano factor, for the symmetric simple 
exclusion process (SSEP). In two dimensions, convergence to the analytical predictions 
was found for large system sizes by assuming a power-law behavior and extrapolating the 
numerical data. In three dimensions no convergence was found. The numerical results 
were, however, obtained for contacts that are not macroscopic in size. Akkermans et 
al. therefore argued that the discrepancy between numerics and theory was caused by 
too small contact sizes with the reservoirs. 

Under certain conditions, the asymptotic current distribution of a one-dimensional 
system that is described by the MFT can be calculated from an additivity principle (AP) 
postulated by Bodineau and Derrida ESI- The validity of this AP has been confirmed in 
several one-dimensional systems, both analytically [25, [261, S7J [28) 29] and numerically 
(2, [2!)][-i0 31|[32]. An interesting question is if one can use the AP to predict the current 
distribution in higher-dimensional systems. This is especially important because many 
experimental systems are higher-dimensional. The results from [24] indicate that it is, 
indeed, possible to do this. So far, only a few studies have addressed this question. Saito 
and Dhar studied heat fluctuations in a deterministic system connected to stochastic 
reservoirs [33] • They found that the AP can predict the current distribution in three 
dimensions, both for diffusive and anomalous heat transport. Hurtado, Perez-Espigares, 
del Pozo, and Garrido confirmed the validity of the AP for the two-dimensional Kipnis- 
Marchioro-Presutti (KMP) model [21 [34] . 

We study numerically the first three cumulants of the current distribution of 
boundary driven generalized exclusion processes (GEPs) [35] . The dynamics is simulated 
using kinetic Monte Carlo (kMC). The simplest case of a GEP is the SSEP, where only 
one particle can occupy each lattice site. In our simulations of the SSEP we consider 
contacts with the reservoirs that are macroscopic in size. Complete convergence of the 
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Fano factor to the analytical prediction of [33] is found in two dimensions. For three 
dimensions the data indicate convergence for large system sizes. We proceed with the 
study of the diffusion coefficient and the current fluctuations in a GEP where maximally 
two (interacting) particles can occupy each lattice site. The first three cumulants of the 
current distribution are calculated by combining the AP with the results from J33]. In 
one and two dimensions the first three cumulants obtained from kMC are in agreement 
with the predicted values. In three dimensions the first two cumulants are in agreement 
with the AP. The statistics for the third cumulant is insufficient for a reliable comparison. 
Because the diffusion coefficient depends on the dimension, the current statistics change 
for different dimensions. The current statistics are independent of the spatial dimension 
for the SSEP and the zero-range process (ZRP). 

The paper is organized as follows. In Section [2] we introduce the quantities that 
are studied. It is explained how to predict the current distribution in any dimension 
from the AP. In Section [3] we present the numerical results for the SSEP. The GEP is 
defined in Section 14.11 The behavior of the diffusion coefficient in different dimensions 
is discussed in Section l4~2l Current fluctuations are studied in Section l4~3l A conclusion 
is presented in Section [5} 


2. Theory 


Consider a one-dimensional system of length L in contact with two particle reservoirs, 
called A and B , at densities P a and P b- The dynamics in the bulk of the system is 
diffusive, i.e., there is no external driving in the bulk. The total number of particles 
that have passed through the system in the time interval [0,f], in the NSS, is denoted 
by Q t • To measure Q t one could, e.g., count the net number of particles entering the 
system from reservoir A. Q t is a stochastic quantity and is described by a probability 
distribution P(Qt)- We study P(Qt ) in the limit t f oo and L f oo. Bodineau and 
Derrida showed that, by postulating an AP, one can calculate the cumulants of P(Qt) 
in a one-dimensional system from the integrals I m [25] [2G] 

rPA 

Im = / D(p)a(p) m - 1 d P . (1) 


' Pb 


( 2 ) 


D(p) is the diffusion coefficient. It is defined by Fick’s first law 

3 = ~D{p)Af, 

where j = (Qt)/t is the average particle flux ((•) denotes the average over P(Q t )), 
and with A p = ps — Pa small enough so that linear response is valid. cr(p ) describes 
equilibrium fluctuations of Q t for large t 




Pa — Pb — P- 


The first three cumulants of P(Qt) are equal to 

(< Qt ) i T 

~r ~ i h: 


(3) 


( 4 ) 
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The ratio of the first two cumulants is called the Fano factor 

F = Urn > ta = !l. 

L^oot^-o o \Qt) i f 

Consider all density profiles pj(x,t'), with 0 < x < L and 0 <t'<t, that lead to 
the same particle flux j. In the long-time limit t f 00 , only the most probable (optimal) 
of these profiles is relevant for the current distribution [53]. The AP is valid as long 
as the optimal profiles are time-independent: pj(x,t') = Pj(x). The point at which the 
optimal profile becomes time-dependent corresponds to a dynamical phase transition 
|22l 1361 EH EE]. For example, for one-dimensional systems on a ring, large fluxes are 
created by traveling waves [36] 137, [38]. One can show from the MFT that a sufficient 
condition on D(p) and cr(p) for the validity of the AP is [55] 

D (p) a "(p) < D\p)a'{p\ Vp. (8) 


Note that (JSJ) is a sufficient but not a necessary condition. 

A qualitative explanation of the AP goes as follows. The system is divided into 
subsystems. Their density profiles are considered to be independent of each other, except 
at the contacts between them. The subsystems should be so small that they are close 
to (local) equilibrium, but yet be large enough to allow for coarse graining. In this case, 
each subsystem has Gaussian current fluctuations around its deterministic behavior 02]) , 
which is completely described by D(p) and cr(p). By calculating the optimal densities at 
the contacts between the subsystems, one finds the cumulant generating function (CGF) 
of the current distribution. From this CGF one can calculate 03}, flH}, OS}- Hence, the 
AP allows one to calculate the current distribution arbitrarily far from equilibrium using 
only the equilibrium quantities D(p) and cr(p). 

We now consider systems in d > 1 dimensions. Fick’s first law is then given by 

J = -D(p)Vp, (9) 


with D(p) a symmetric d x d matrix. If the diffusion is isotropic, which is the case 
considered here, one can write D(p) = D d (p)I d , with D d (p ) a scalar function depending 
on the dimension. A sufficient condition that excludes the possibility of a dynamical 
phase transition is OS]) with the scalar functions D d (p) and a d (p) [22]. 

Akkermans et al. studied current fluctuations in higher-dimensional diffusive 
systems The shape of the system and the contacts with the reservoirs are taken 
arbitrary, but macroscopic in size. If the optimal density and current profiles are time- 
independent, the MFT predicts that the CGF of the system in d dimensions /id(A) equals 

m 




( 10 ) 
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with k a constant that depends on the shape of the system and shape of the contacts 
with the reservoirs. The calculation of k is explained in Appendix D pflX) is the CGF 
of a one-dimensional system described by Dd(p) and <Jd(p)- Since one assumes that the 
optimal density and current profiles are time-independent, pflX) can be calculated from 
the AP, by using Dflp) and (Jd{p) in (JH) • 


3. Symmetric simple exclusion process 

The SSEP is a stochastic lattice gas where particles interact by exclusion, i.e., each 
site can contain maximally one particle. Each particle attempts to hop to its nearest 
neighbors with unit rate. A hopping attempt is successful if the site is empty. The 
distance between two sites is equal to one. For the SSEP D(p) = 1 and a(p) = 2p(l — p) 
in any dimension. (JHJ) is therefore always satisfied. We consider reservoirs with densities 
Pa = 1 and ps = 0. A calculation from the AP [25] or an exact microscopic derivation 
[ 39] shows that F = 1/3 in one dimension. Since D(p) and er(p) are independent of 
the dimension, F = 1/3 in any dimension. It is, however, important that the size of 
the contacts scales with the system size, thereby maintaining a finite fraction of the 
boundary in contact with the reservoirs. The numerical computation of the Fano factor 
in [24] was performed for systems where this scaling is absent. We present simulations 
in which the contacts do scale with the system size. 


The dynamics is simulated using a kMC algorithm, cf. Appendix A How the 


Fano factor is computed from the simulation data is explained in |Appendix B| In two 
dimensions we consider squares of size L x L and in three dimensions cubes of size 
L x L x L. The contact between the system and the reservoirs is modeled as lattice 
sites whose densities are fixed and uncorrelated from the rest of the system, as in [ 24] . 
The shape of the contacts is illustrated in Figured! 

The numerical results for the Fano factor are presented in Figure [2aJ For two 
dimensions the Fano factor has converged to 1/3 at L « 40. This extends the 
extrapolation presented in Figure 3 of [24]. We determined numerically that k ~ 0.663L 


for the geometry in Figure [2a] cf. Appendix D The average current indeed converges to 
L(Qt)/t ~ 0.663L, compared to L(Q t )/t = 1 in one dimension (data not shown). For 
three dimensions convergence to 1/3 is not yet attained at L — 15. However, the data 
indicate convergence to 1/3 for larger system sizes. For the same distance L between the 
two reservoirs, the Fano factor in three dimensions is lower than in two dimensions. One 
therefore expects convergence before L — 40 in three dimensions. In Figure [2b] we plot 
the Fano factor in three dimensions as a function of 1/L 2 . There is no specific reason 
to assume that this is the correct convergence law. We choose this scaling because we 
want to compare our results with Figure 4 of [24]. A 1/L 2 fit indicates an L —> oo limit 
of F = 0.3344, with one-sigma error bar a = 0.0018. The fit was performed using the 
method of least squares with weighted error bars in. The extrapolation is in agreement 
with the expected value of F — 1/3. Our numerical results validate the conjecture from 
m that the observed discrepancy between numerics and theory is caused by too small 
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(a) 
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Figure 1: The type of contacts used for the SSEP in two dimensions (a) and in three 
dimensions ( b ). The black dots are sites with a particle density of 1 (A) or 0 (B), whose 
state is uncorrelated from the rest of the system. In two dimensions, 1/2 of the lower 
left is connected to reservoir A and 1/2 of the upper right is connected to reservoir B. 
In three dimensions, 2/3 of the lower left is connected to reservoir A and 2/3 of the 
upper right is connected to reservoir B. 




(a) 


(b) 


Figure 2: (a) The Fano factor with one-sigma error bars for the SSEP, for squares Lx L 
and cubes L x L x L as depicted in Figure Q] The lines are a guide to the eye. The 
two-dimensional results show a convergence to 1/3 at L ~ 40. The three-dimensional 
results have not yet converged. (6) The three-dimensional data as a function of 1/L 2 for 
L > 9. The thin black line is a 1/L 2 fit using the method of least squares with weighted 
error bars. The thick black lines are one-sigma error bars on the L —y oo limit predicted 
by the fit. 


contacts with the reservoirs. 
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4-1. The model 

Recently, we have studied the diffusive behavior of a lattice model of interacting particles 
mmm- The original motivation was the study of diffusion in nanoporous materials 
[44] . Some of these materials have a structure consisting of cavities connected by narrow 
windows [45], as illustrated in one dimension in Figure [3] Each cavity can be identified 
as a lattice site and can contain between 0 and n max particles. The distance between 
two lattice sites is taken equal to one. The length of the system is the distance between 
the two reservoirs L = N + 1, with N the number of cavities. A cavity containing 
n particles has an equilibrium free energy F(n ) that depends solely on the number of 
particles it contains. If the system is in equilibrium at chemical potential /i and inverse 
temperature fl = (/c^T) -1 (with kb the Boltzmann constant), the probability to observe 
n particles in any cavity is equal to 

pM =[Z(p)}- 1 (11) 

with Z the grand-canonical partition function 

Tim ax 

z{/x) = j2 e ~ mn) ~ >in] - ( 12 ) 

77.—0 

Averages over the equilibrium distribution (fill) are denoted by (•), e.g., 

TT-m ax 

( 13 ) 

77=0 

(Whether (•) denotes the average over p e fl(p) or P(Qt) is always clear from the context.) 
Particles jump from a cavity containing n particles to a cavity containing m particles 
with rate k nm . These rates obey local detailed balance p e flp e flk nrn = p^ +1 p® q _ 1 fc m+ i !n _i. 
The reservoirs are modeled as cavities whose probability distribution is uncorrelated 
from the rest of the system. The rates at which a reservoir cavity at chemical potential 
H adds (/c+) or removes (k~) one particle from a cavity containing n particles are 

Tim ax Timax 1 

k n = ^ ^ ^mnPm(h); k n = ^ ^ k nm p r ^(p). (14) 

777=1 777=0 

This model is a GEP [35] with a stochastic thermodynamical interpretation for the 
equilibrium statistics and dynamics. When defined like this it is an adequate model for 
the understanding of the equilibrium and diffusive behavior of particles in nanoporous 
materials [4111421 143] . For n max = 1 the model reduces to the SSEP. A zero-range process 
[46 ] is defined by rates that only depend on the cavity from which the particle jumps. 
Hence, one finds a ZRP for n max = oo and k nm = k n . 

In the following, we fix the parameters to n max = 2 and fl = 1. The free energy 
can be written as F{n ) = In n! + cn + /(n), with c a constant j jl ■ 02]. The first 
term accounts for the indistinguishability of the particles. The linear term cn is the 
ideal gas contribution. f(n) is nonzero because of particle interactions, and is called 
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Figure 3: A lattice model of a nanoporous material. Each cavity (upper drawing) is 
mapped to a lattice site (lower drawing) and contains between 0 and n max particles (here 
n ma x = 2). On the boundaries the system is connected to cavities that are uncorrelated 
from the system. A cavity containing n particles has equilibrium free energy F[n). 


the interaction free energy. Note that a linear term in F(n) is equivalent to adding a 
constant to the chemical potential p (1111) . A linear term does therefore not influence 
the equilibrium statistics at a given particle concentration. The rates we consider are 

k nm — ne -[/( n - 1 )+/( m + 1 )-/( ri )-/( m )]/ 2 _ (]_ 5 ) 


It is clear that a linear term in F{n) (or f(n)) also does not influence these rates. Hence, 
we can rescale F{n) so that /(0) = /(1) = 0 without loss of generality. All possible 
interactions are then described by /(2). In the following we consider attractive particles 
/(2) < 0. In this case, the form of the transition rates (1T51) can be rationalized from 
transition-state theory [32]. Furthermore, for this choice of rates the diffusive behavior 
agrees with experiments of attractive particles in nanoporous materials m- 

For an isothermal system, which we consider here, D(p ) and cr(p) are related by 
the following fluctuation-dissipation relation [47] 

a(p) = 2 k b Tp 2 K(p)D(p), (16) 


with n(p) the isothermal compressibility. One knows from statistical physics that n(p) 
can be written as 


K(p) = fl 


V (n 2 ) - (n) 2 
(n) (n) 


(17) 


with V the volume in which the average ( n) and particle fluctuations ( n 2 ) — ( n ) 2 are 
measured. Because particles in different cavities do not interact, one can take the 
averages over one cavity, V = 1 and p = (n). One then finds for cr(p) (fT6l) 


a(p) = 2(( n 2 } - (nf)D(p). 


(18) 


Regarding notation, since p = (n) we use p and (n) interchangeably. Also, averages 
(•) are a function of the chemical potential of the reservoirs. These can, however, be 
straightforwardly converted to densities via ([T3lh In this paper we write everything as 
a function of the density. 
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From (fTgp one finds that I m ([I]) can be written as 

[■Oa 


I'm. 


D((n)r [2((rt 2 } - <n) 2 )] m_1 d(n), 


Wi 


(19) 


where (u)a and ( n)s are the average number of particles in, respectively, reservoir cavity 
A and B. One can compute I m by numerically simulating D((n)) and analytically 
calculating ( n 2 ) — ( n ) 2 from p e fl(p). 


f.2. Diffusion coefficient 

We have studied D(p ) in this model both numerically and analytically [4Tj|42, 43]. From 
these studies one can conclude that D(p) is, in general, influenced by correlations (see 
also BE!). Since the effect of correlations changes and is actually expected to diminish 
with increasing dimension, the function D(p) depends on the dimension | [42l [48] . If the 
effect of correlations upon the diffusion is completely neglected one can show that D(p) 
is given bv [ill 142] 

(k) 


D{p) — 




(n 2 ) — (n) 5 


( 20 ) 


This result is valid for a (hyper)cubic lattice in any dimension. Because one arrives 
at (l20|i by neglecting all correlations, it could be argued that in the limit of infinite 
dimension D(p) converges to (120T) . Although we do not have a rigorous proof of this 
statement, it is confirmed by numerical evidence given below (see also 03). We therefore 
denote the results that are calculated from (120|) as the d —* 00 limit. Note that in this 
limit the integral (TT9l) can be calculated analytically. 

The uncorrelated result (T20|) is exact for the SSEP (n max = 1), which is easily 
checked by using that pff = p and pff = 1 — p. It is also the same in any dimension j49j. 
(|20l) is also exact for the one-dimensional ZRP [43]. Since the particle distribution in 
the NSS factorizes in any dimension for the ZRP [50], the calculation from [43] can be 
straightforwardly extended to higher dimensions to show that D(p) is independent of the 
dimension. To our knowledge, these are the only two cases where the uncorrelated result 
is exact for GEPs. It is, then, no surprise that D(p) is independent of the dimension. 

We consider now /(2) = —2.5. This is a concave /(n), signifying attractive particles 
[141] . We choose this interaction because correlations have a large effect for attractive 
particles. In Figure [4] we plot D(p) in one, two, three, and infinite dimensions. We 
refer to Appendix C for details on the simulations. D(p) appears to converge with 
increasing dimension towards the d —> 00 result (j20lh The diffusion coefficient as a 
function of the dimension for ( n ) ~ 0.51 and (n) ~ 1.49 is shown in, respectively, 
Figures [5k and [5fc>. The behavior is well approximated by a 1/d dependence. Figure 
[5b shows the same quantity for the interaction /(2) = 0 at ( n) = 1 (data from [48]). 
Also here an approximate 1/d dependence is found. This dependence can be understood 
as follows. Correlations are the result of memory effects in the environment [42]. The 
strongest contribution comes from the increased probability that a particle jumps back 
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Figure 4: D(p) for /(2) = —2.5 and n max = 2 in one, two, three, and infinite dimensions. 
The error bars are smaller than the symbol sizes. 



Figure 5: The diffusion coefficient as a function of the dimension, for n max = 2. The 
data are normalized w.r.t. the analytical uncorrelated result fl20]h which is denoted by 
.D(oo). The black circles are from kMC simulations and the red squares are (1201) . The 
error bars are smaller than the symbol sizes. 1/d fits were performed with the method 
of least squares. In all three cases this fit provides a good estimate for the diffusion 
coefficient at infinite dimension, with a relative error (.Dfi t (oo)/-D(oo) — 1) of a) 0.3 %, 
b) 2.0 %, and c) 0.07 %. 


to its previous position. The probability to do so is approximately 1/2 d as there are 2d 
neighboring cavities. This simple argument indeed suggests that the effect of correlations 
will decrease approximately as oc 1/d. 
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Figure 6: Plot of D(p), a(p), and D'(p)a'(p) — D(p)a"(p) for /(2) = —2.5 and 
w m ax = 2 in the limit d —* oo (120T) . The sufficiency condition (JHJ) is satisfied if 
D'{p)(j'{p) — D(p)a"(p) > 0 for all p. 


4-3. Current fluctuations 


The sufficiency condition (J8]) is not satisfied for /(2) = —2.5, as shown in Figure |6] for 
d — > oo (l20j) . The numerically simulated H(p)’s do not give smooth results for ([8]), since 
one has to calculate the second derivative of an interpolated function. The qualitative 
behavior of (JHJ) for finite dimensions is, however, the same as for d —> oo. Starting from 
([20]) . one sees that (J8J) does not hold for many GEPs. One can show analytically that all 
GEPs with n max = 2 and /(2) <0 do not satisfy (J8l). Numerically, one folds that GEPs 
with n max = 2 and /(2) > 2.917 also do not satisfy (jHJ). Although (JHJ) is not satisfied for 
the parameters considered here, we expect that the AP is still valid. Dynamical phase 
transitions have only been observed for closed systems [22, 36| [3?, [38], not boundary 
driven ones [30], [31, [32j. Also, dynamical phase transitions do not occur for currents 
close to the average current [23]. Currents created by time-dependent density profiles, 
if any, are therefore highly unlikely, and their influence on the first three moments of 
the current distribution is expected to be negligible. 

We study the current statistics for /(2) = —2.5 and reservoir densities ( h)a = 
Umax = 2 (pa = oo) and (u)b = 0 (ps = —oo). Plots of L(Q t )/t = I\ and 


L{Qt)c/t = h/Ii as a function of the length are shown in, respectively, Figures I7al 
and m The values predicted by the AP are given by lines, which are the one-sigma 
error bars. These error bars arise from the error bars on the simulated D(p)’s. Values 
from direct numerical simulations are also given with one-sigma error bars, as explained 


in 


Appendix B 


Let us first consider the one-dimensional data. We estimate convergence in length 


at L pc 175. How we check for convergence in time is explained in Appendix B The 
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N = L — 1 
(a) 


N = L-l 
(b) 


Figure 7: (a) L{Q t )/t = h and (b) L(Q 2 ) c /t = I 2 /h, for /(2) = -2.5, n max = 2, and 
different lengths in one, two, and three dimensions. The lines are predictions from the 
AP, and represent one-sigma error bars. The points with one-sigma error bars are from 
a direct simulation of the current. In two (three) dimensions, the directly measured 


cumulants are divided by L y (L y L z ), see Appendix D 


value for I \, cf. Figure [HaJ is taken from the highest considered length in Figure [Taj 
To achieve a good statistics for the second and third cumulant, we have performed an 
extensive simulation for length L = 251. The simulated values for L(Q 2 ) c /t = h/h, 
F = h/df, and L(Q\) c /t for this length are given in, respectively, Figures l8bl l8cl and [Ml 
I\ from the AP is slightly higher than the directly simulated value (/ Ap /If m fa 1.0018). 
The most likely reason for this is that the simulated D(p) slightly overestimates the 
real D(p). The diffusion coefficient should be measured in the limit of an infinitely 
small concentration gradient, while of course the simulations are performed at a finite 
concentration gradient. Similarly, one should in principle simulate an infinitely long 
system, so that all boundary effects have disappeared. Both approximations cause 
the numerically simulated D(p) to overestimate the real value J-JEJ. Furthermore, to 
calculate J AP one has to interpolate the simulated points of D(p), and then integrate 
this interpolated function. This could introduce a small numerical imprecision. Since 
the relative difference is less than 0.2% we consider this result a very good agreement 
between / AP and /f m . Also the variance and the Fano factor are in very good agreement 
with the value from the AP: (J 2 AP /J 1 AP )/(/ 2 sim // 1 sim ) « 1.0007 and F AP /F sim « 0.9989. 

Figure [870 shows the third cumulant. Although the error bars are significantly larger 
compared to the first two cumulants, the data indicate agreement between the AP and 
the directly simulated values. Finally, we plot P(Qt) obtained from kMC together with 
the Gaussian prediction from the first two moments of the AP in Figure [9J The small 
error on I\ from the AP is noticeable for determining (Qt)- When using the simulated 
(Qt), one sees that P(Qt) is well approximated by a Gaussian. Indeed, the skewness of 
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d 



d 


(a) 


(b) 



d d 


(c) 


(d) 


Figure 8: (a) h, (b) h/h, (c) F, (d) L(Qf) c /t for /(2) = -2.5 and n max = 2 as a 
function of the dimension. Predictions from the AP are denoted by blue error bars 
without symbol. The limiting case d —» oo is shown as a black line. Direct numerical 
simulations are denoted by black diamonds. In two (three) dimensions, the directly 
measured cumulants are divided by L y (L y L z ), see Appendix D Note that the error bar 


at d = 3 in (d) spans approximately 7 times the whole y axis. 


P(Qt ) is small (Ql) J (Q 2 )^ 2 ~ 0.034, i.e., P(Q t ) is almost symmetric. Although the 
difference is small, one observes that for Q t < 320 the simulated P(Qt) is consistently 
lower than the Gaussian, while for Q t > 365 it is consistently higher. 

We now discuss the higher-dimensional systems. In contrast to the SSEP, all sites 
at the boundaries are in contact with the reservoirs. If periodic boundary conditions 
are imposed in the y direction, D(p) converges in two dimensions to the L y f oo limit 
at L y k, 3. In the simulations we take L y = L z = 5 with periodic boundary conditions. 
The diffusion coefficient is simulated for the same concentration gradients and length in 
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Qt 



Qt 


(a) (b) 

Figure 9: P(Qt) from kMC (red) in one dimension for n max = 2, /(2) = —2.5, L = 251, 
and t = 8.10 4 . A Gaussian distribution (black) with average and variance predicted by 
the AP (a) and the simulated average and variance from the AP (b) is also plotted. The 
data is well approximated by a Gaussian distribution. 


the x direction as for the one-dimensional case. The different coupling to the reservoirs 
compared to the SSEP is done for numerical reasons. The program for the GEP is too 
slow to simulate a convergence in both the x direction and y (z) direction. The periodic 
boundary conditions employed here are equivalent to the L y (L z ) f oo limit. All sites 
at the boundaries are coupled to the reservoirs because this gives the highest particle 
flux. The higher the particle flux, the better the current statistics for a given simulation 
time. If all boundary sites are connected to the reservoirs k = L y and k = L y L z 


in, 


respectively, two and three dimensions. This is explained in Appendix D 


For two dimensions we assume convergence in length at L ~ 120. The error on I\ 
is comparable to the one-dimensional case (/^ p //f m ~ 1.0010). The second and third 
cumulants are determined from extensive simulations at length L = 121. The variance 
and Fano factor are slightly underestimated by the AP: (/^ LP //j 4p )/(/| im //f m ) « 0.9982 
and F AP /F sim « 0.9971. We consider this a very good agreement between the direct 
simulations and predictions from the AP. The relative difference is less than 0.3 %, and 
all quantities show a large overlap within their error bars. The third cumulant is also 
compatible with the AP prediction, although the error bar on the directly simulated 
value is rather large. The shape of P{Qt ) is similar to the one-dimensional case (data 
not shown). 

For three dimensions the simulation times become much longer. We therefore only 
simulate the current for systems of length L = 101 and L = 121. Since the two- 
dimensional system has converged at L = 121, one can safely assume that this is also 
the case for the three-dimensional system. The cumulants from Figure [8] are calculated 
for L = 121. The average, variance, and Fano factor are correctly predicted by the AP. 
















Current fluctuations in boundary driven diffusive systems 


15 


There is insufficient data to achieve a reliable estimate for the third cumulant. The 
obtained value shown in Figure [8d] agrees well with the AP, but the error bar is very 
large: L{Q\) c /t + a = 0.438 and L(Qf) c /t — a = —0.110. The shape of P(Qt) is similar 
to the one-dimensional case (data not shown). 

5. Conclusion 

To conclude, we have studied numerically current fluctuations in the symmetric simple 
exclusion process (SSEP) and a generalized exclusion process (GEP). For the SSEP we 
find that the Fano factor is independent of the spatial dimension and (macroscopic) 
shape of the contacts with the reservoirs. For the GEP our numerical simulations are in 
agreement with the predictions from the AP combined with the MFT [23]. In one and 
two dimensions agreement is found for the first three cumulants. In three dimensions 
the first two cumulants agree with the AP, while the statistics for the third cumulant 
are insufficient for a reliable comparison. The diffusion coefficient, and as a result the 
current statistics, depends on the dimension. Only for the SSEP and the ZRP is the 
diffusion coefficient independent of the dimension. 

A more precise numerical determination of the diffusion coefficient from Fick’s first 
law is computationally very time consuming, at least using the methods presented here. 
It would therefore be of interest to find exact analytical results for the diffusion coefficient 
for the GEP. Another interesting question concerns the simulation of higher moments 
of the current distribution. This could be achieved using a sophisticated Monte Carlo 
algorithm to simulate rare events, see e.g. [51[52, 53]. Both the SSEP and the ZRP 
satisfy the sufficiency condition for the validity of the AP (J8]). However, many GEPs 
do not satisfy (J8]). Hence, one might observe deviations from the predictions of the AP 
for large current fluctuations. An analysis of the optimal density profiles, before and 
(possibly) after the dynamical phase transition, is also of interest. 

The quantities D(p) [53] and a(p) [55] are experimentally accessible in nanoporous 
materials. The average particle flux through a system in contact with two particle 
reservoirs can also be measured [56]. If it is possible to measure the variance of the 
particle flux with a good precision, these techniques present an opportunity for an 
experimental verification of the additivity principle and, therefore, the macroscopic 
fluctuation theory. 
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Because we consider pA — 1 and ps = 0 for the SSEP, all transition rates are equal to one 
(also at the boundaries). All n possible transitions are stored in a list. A random integer 
between 0 and n — 1 decides which transition takes place. The time between two events 
is taken from the distribution p(t) = nexp(— nt). For the GEP with n max = 2 there are 
12 different rates (four in the system and four at the contact with each reservoir). Since 
this is a small number, one can use the algorithm described by Schulze EZt For a fixed 
number of Monte Carlo steps, the computation time of both algorithms is constant for 
different system sizes. 


Appendix B. Data analysis 


The current fluctuations are measured as follows. First the system is allowed to relax 
to its steady state, after which we put the time at 0. The net number of particles that 
have entered the system between time 0 and t is denoted by Qt,i- The net number of 
particles that have entered between time t and 2 1 is denoted by Qtp, and so on. In 
the simulations Qt is determined by measuring the particle current at the left and right 
boundary. One then has a list {Qt} with A) elements. The average is equal to 
N t 

Ol=Y.Qc/ N i- (B- 1 ) 

i— 1 


For large A) the average Q t is a good approximation for the average (Qt) over P(Q t ). 
The sample variance is equal to 


Ni 

St 2 = £(<?M-® 2 /W-l)- (B.2) 

2—1 

For large A), S‘f converges to (Qf) — (Qt) 2 - 

The one-sigma error bar on Q t is equal to (assuming the Q t f s are independent 
identically distributed variables) 


<r= f S?/N,. 

The variance of Sf is equal to 


Var(*S' t 2 ) 



(B.3) 


(B.4) 


with <74 = ((Qt — (Qt)) 4 ) the fourth central moment of P(Qt ) (see for example exercise 
7.45 in [SS])- We estimate o by (IB .31) . We do not estimate cr 4 directly from the simulation 
data, because our data do not allow for an accurate prediction of the fourth moment. 
Rather, we use the prediction for cr 4 from the AP [25] . One-sigma error bars on S\ are 
equal to [Var (S) 2 )] 1//2 . Except for the third cumulant, all other error bars are obtained 
from addition and multiplication of Q t and Sf. The rules for finding these error bars 
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can be found in e.g. |2E|. The Fano factor is calculated by F{t) = Sf /Q t . The error bar 
on the third cumulant is found by bootstrapping the simulated data. 

By adding the currents pairwise Q t ,i + Qt,i +1 (with i odd), one can calculate Q2t 
and S% t for the time interval 2 1 (with Nfl 2 points), and so on. We study the Fano factor 
F(nt) for 1 < n < 6. 

We now explain how we check if the data have converged in time. For clarity we 
consider the specific example of the two-dimensional SSEP at L = 40 with t = 2.10 4 . 
The autocorrelation (AC) of Q t)i and Qt,i+i is 




AC = 


XA (Qt,i — Qt 


(B.5) 


The AC is plotted in Figure IBlal together with the critical values (CVs) to reject the 
null hypothesis that AC = 0 at 95 % significance level. All points are smaller than the 
CVs. The point at n — 1 is, however, very close to the lower CV. This suggests that 
there is still a non-negligible AC for times It. Indeed, for small times the AC is always 
negative. For large times, when the Q t fl s are uncorrelated, the AC fluctuates close to 
zero. The scale of “close to zero” is determined by the CVs. 

The Fano factor F{nt) is plotted in Figure IBlbl F(lt) is slightly higher than the 
other 5 points, indicating again that there is not yet convergence in time. The first 
two point that are converged in time are F(2t) and F{3t). A plot as a function of the 
number of simulated points IVj for F(3t) is shown in Figure lB2l After N{ ~ 25.10 4 the 
data fluctuate around the end value F fina i, indicating a good convergence for F(3t). The 
average of F(2t) and F(3t) is taken as the final data point (as plotted in Figure l2alh 
For most points, the first two converged values are averaged to calculate the final result. 
If computation times are exceedingly long, such as for the SSEP in two dimensions for 
L = 50, only the first converged point is taken. In this case this point is F{2t). F(3t) 
has not yet converged as can be seen from a graph similar to Figure lB2l This explains 
the large error bar for L = 50 compared to the other points for the two-dimensional 
SSEP. 


Appendix C. Simulation of diffusion coefficient 

D(p) is simulated for 30 concentrations, see (41, :42] for details on the simulations and 
calculation of the error bars. In this paper the length in the x direction is L = A+l = 16 
in two and three dimensions. In one dimension the analysis was performed for L = 21 
and L = 16. The predicted values of I\ were the same up to a relative difference of 
0.006%. The data in the paper are for L — 21 in one dimension. The concentration 
gradient for low and high concentrations is taken between A p = 0.05 and A p = 0.03. For 
the other concentrations we take A p = 0.06. The values at p = 0 and p = n max can be 
calculated analytically: 71(0) = 1 and D( 2) = 2. An approximation for the continuous 
function D(p) is achieved by interpolating these 32 points (using the “Interpolation” 
function of Mathematica). For concentrations smaller than p ps 0.04 and higher than 
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Figure Bl: The two-dimensional SSEP with t = 2.10 4 , L = 40, and the geometry of 
Figure [Ta] (a) (circles) Autocorrelation (IB. 51) . (red squares) Critical values to reject the 
null hypothesis AC = 0 at 95 % significance level, (b) (circles) F{nt). (dashed line) 
Average of F(2t) and F(3t). This is the value of the data point in Figure l2al 



Figure B2: The two-dimensional SSEP with t = 2.10 4 , L = 40, and the geometry of 
Figure[lal (thick black line) F(3t) after iVj simulated points, (thin grey lines) one-sigma 
error bars, (dashed line) final value of F(3t). 


p 1.96 the interpolated values are higher than the uncorrelated result ([20]) . Since 
we know that correlations lower D(p), we consider the uncorrelated results for these 
concentrations instead of the interpolated function. 
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The CGF /id(A) of a d-dimensional system is equal to (cf. the last equation in [23] ) 


PdW — 


' d— 2 


dr 


x [UMA)]. 


(D-l) 


Hi( A) is the CGF of a one-dimensional system described by D d (p ) and a d (p )• Consider a 
rectangular system of length L x and height L y . All sites at x = 0 are coupled to reservoir 
A and all sites at x = L x are coupled to reservoir B. L is the typical domain size, which 
we take equal to L x . v(x, y) is a function on the domain 0 < x < 1, 0 < y < L y /L x , that 
satisfies the Laplace equation A v(r) = 0, with v(0,y) = 0, v(l,y) = 1, and Neumann 
boundary conditions otherwise. For the geometry we consider it is straightforward to 
show that v(x, y) = x. One then finds 

pLy/L X 

dxdy x [L x yi( A)] = L y pfl\). (D.2) 

o Jo 


/12(A) — 



The calculation for the same geometry in three dimensions shows that /x 3 (A) = 
LyL z pi (A). 

The density p(x, y) can be found from the one-dimensional profile p\ (x) (equation 
(33) in [24]) 

p(x, y ) = Pi(v(x, y)) = pflx). (D.3) 


Note that the only assumption required for these results is the time-independence of 
the optimal density and current profiles. In the study of the two-dimensional KMP 
model with all the boundary sites connected to reservoirs |2, [33], one made the extra 
assumption that the optimal current profile is constant jj(r) = J. This extra assumption 
is unnecessary: it can be derived from the time-independence of the optimal profiles and 
the MFT. Indeed, in one dimension time-independent profiles imply a constant current 
profile. A constant current profile in two dimensions follows from (ID. 3D . Note that for 
more general couplings to the reservoirs, such as in Figure [U the optimal current profile 
need not be constant. 

We have solved numerically the Laplace equation for v(r) for the domain in Figure 
l2al One finds / 12 (A) ~ 0.663L/ii(A). This agrees with our kMC results, as discussed in 
Section [3] 
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